home *** CD-ROM | disk | FTP | other *** search
/ CD Concept 6 / CD Concept 06.iso / mac / UTILITAIRE / RLaB / testmatrix / dual.r < prev    next >
Text File  |  1994-12-27  |  2KB  |  76 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Dual vector with respect to Holder p-norm.
  4.  
  5. // Syntax:    dual ( X , p )
  6.  
  7. // Description:
  8.  
  9. //    dual(X, p), where 1 <= p <= inf, is a vector of unit q-norm
  10. //    that is dual to X with respect to the p-norm, that is, norm(Y,
  11. //    q) = 1 where 1/p + 1/q = 1 and there is equality in the Holder
  12. //    inequality: X'*Y = norm(X, p)norm(Y, q).
  13.  
  14. //    Special case: DUAL(X), where X >= 1 is a scalar, returns Y
  15. //    such that 1/X + 1/Y = 1.
  16.  
  17. //      Called by PNORM.
  18.  
  19. //    This file is a translation of dual.m from version 2.0 of 
  20. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  21. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  22.  
  23. //-------------------------------------------------------------------//
  24.  
  25. dual = function ( x , p )
  26. {
  27.   local(p)
  28.   
  29.   if (max (size (x)) == 1 && nargs) 
  30.   { 
  31.     p = x;
  32.   }
  33.  
  34.   // The following test avoids a `division by zero message' when p = 1.
  35.   
  36.   if (p == 1)
  37.   {
  38.     q = inf();
  39.   else
  40.     q = 1/(1-1/p);
  41.   }
  42.  
  43.   if (max (size (x)) == 1 && nargs == 1)
  44.   {
  45.     y = q;
  46.     return y;
  47.   }
  48.  
  49.   if (norm (x,"i") == 0)
  50.   {
  51.     return x; 
  52.   }
  53.  
  54.   if (p == 1)
  55.   {
  56.     y = sign (x) + (x == 0);     // y(i) = +1 or -1 (if x(i) real).
  57.  
  58.   else if (p == inf()) {
  59.  
  60.     xmax = max (abs (x));
  61.     f = find (abs (x) == xmax);
  62.     k = f[1];
  63.     y = zeros (size (x));
  64.     y[k] = sign(x[k]);           // y is a multiple of unit vector e_k.
  65.  
  66.   else  // 1 < p < inf.  Dual is unique in this case.
  67.  
  68.     x = x/norm (x,"i");          // This scaling helps to avoid under/over-flow.
  69.     y = abs(x).^(p-1) .* ( sign(x) + (x == 0) );
  70.     y = y / norm (y, q);         // Normalize to unit q-norm.
  71.     
  72.   } }
  73.  
  74.   return y;
  75. };
  76.